Типы пространственных данных;
Особенности геоданных;
R и пространственный анализ;
Практика - находим геоданные разных форматов и учимся с ними работать;
Что можно почитать и посмотреть по теме.
а. Векторные данные
Геометрии состоят из точек. Конвенционально x = широта, y = долгота. Кроме этого могут использоваться z = высота, M = мера изменчивости точки;
Виды геометрий:
+++ GEOMETRYCOLLECTION +++
Форматы записи данных:
Well-Known Text (WKT) - Запись точек вектором из геометрий
Well-Known Binary (WKB) - Запись координат в бинарных значениях. Используется в базах данных (т.к. увеличивает скорость работы с данными), но нечитаем и непонятен для человека.
Посмотрим в первом приближении:
library(sf)
library(rnaturalearth)
world <- ne_countries(returnclass = "sf")
world <- world[1:5,c(4, ncol(world))]
# Для иллюстрации
world_1 <- world
world_1$geometry <- as.character(world_1$geometry)
as.data.table(world_1)
b. Растровые данные
Матрица значений пикселей пространственной области. В растрах хронят космические снимки и базы различных геологических данных. Также как и вектор - предназначены для определённого круга задач;
Основной формат: .tif
Основные библиотеки для обработки в R: raster и stars
Сравним:
На этом и последующих занятиях мы будем говорить только о векторных данных.
Главное, нужно помнить, что векторные пространственные данные - это (чаще всего) просто точки на плоскости! Создадим собственные пространственные данные:
p1 = st_point(c(7,52)) # Рисуем точку 1
p2 = st_point(c(-30,20)) # Рисуем точку 2
sp_obj = st_sfc(p1, p2, crs = 4326) # Преобразуем в пространственные данные
plot(sp_obj)# Вы восхитительны!…
Вернёмся к векторным геометриям:
world
Coordinate Reference Systems (CRS) / Пространственная привязка и её проекции
Что это такое?
Для иллюстрации обратимся к прекрасной работе Тасс.
Картинки чтобы окончательно закрепить идею проекции.
Как посмотреть проекцию в CRS?
st_crs(world)## Coordinate Reference System:
## User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## wkt:
## GEOGCS["WGS 84",
## DATUM["WGS_1984",
## SPHEROID["WGS 84",6378137,298.257223563,
## AUTHORITY["EPSG","7030"]],
## AUTHORITY["EPSG","6326"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4326"]]
На самом деле всё не так сложно.
Пространственный объект можно перевести в другую CRS сразу и быстро, задав через st_transform() нужную проекцию одним из четырёх форматов (просто гуглите то, что нужно для вашего пространства):
Например:
Afg = world[1,]
plot(Afg)# Change crs
Afg <- st_transform(Afg, 2264)
plot(Afg)# Change crs
Afg <- st_transform(Afg, "+proj=longlat +datum=WGS84 +no_defs")
plot(Afg)В геоданных есть ещё множество подводных камней, но этот - ключевой и его необходимо запомнить.
Почему R?
Библиотеки для работы с геоданными в R
(несть им числа, но основные)
Немного истории:
Что мы будем делать - делать карту мороженного в Израиле! :)
Я исхожу из того, что вы знаете dplyr или data.table на базовом уровне, но если у вас вызывает затруднения именно манипуляции с данными - обязательно задавайте вопросы!
Мы отработаем работу с данными на трёх “классических” форматах для векторных данных: csv, geojson, shape-file.
A. Датасет с мороженками в csv
Загрузим нужные библиотеки:
library(sf)
library(geojsonsf)
library(geojsonio)
library(readr)
library(data.table)Загрузим данные. К сожалению, в R есть известная проблема с кодировками и чтобы долго не решать проблему с крокозябрами - уберём колонки названий точек мороженного на иврите:
isr_icecream <- fread('https://raw.githubusercontent.com/AmitLevinson/Datasets/master/golda/golda_locations.csv', encoding = 'UTF-8')
isr_icecream[,1:2 := NULL]
isr_icecreamПревратим в sf тип данных:
isr_icecream <- st_as_sf(isr_icecream, dim = "XY", remove = T, na.fail = F, coords = c("lon", "lat"), crs = "+proj=longlat +datum=WGS84 +no_defs")
isr_icecreamplot(isr_icecream)
B. Границы Израиля в Shape-file
Идём сюда, ищем Израиль и загружаем архив. Что мы видим внутри?
.shp – главноеый файл с геометриями
.dbf
.shx
.prj
(но все они важны!)
По этой причине, когда вам нужно переслать shape-file отправляйте архив из всех четырёх файлов.
Загрузим его и отобразим:
isr_border <- st_read('/cloud/project/ISR_adm/ISR_adm0.shp')## Reading layer `ISR_adm0' from data source `/cloud/project/ISR_adm/ISR_adm0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 34.26801 ymin: 29.49708 xmax: 35.90094 ymax: 33.36403
## CRS: 4326
isr_borderplot(st_geometry(isr_border))
C. Точки городов Израиля в GeoJson
Json с геоданными. Достаточно тяжёлый с точки зрения размещения данных, но часто встречающийся формат.
И здесь стоит отвлечься и ознакомиться с “Википедией”" от мира пространственных данных - Open Street Map.
Если это “Википедия”, то с неё можно выкачать данные? Да. По тэгам…
…
Тэги городов.
Скачаем GeoJson с городами. Загрузим его. Внутри будет огромное количество колонок, которые сейчас нам не нужны. Оставим первую и последнюю:
isr_city <- st_read('/cloud/project/city.geojson')## Reading layer `city' from data source `/cloud/project/city.geojson' using driver `GeoJSON'
## Simple feature collection with 51 features and 225 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 33.80459 ymin: 30.61197 xmax: 36.56875 ymax: 32.92817
## CRS: 4326
isr_city <- isr_city[,c(1,ncol(isr_city))]
isr_cityplot(isr_city)
Альтернативный способ через пакет osmdata (спасибо Philipp Pyshny за наводку)
library(osmdata)
admin_osm = opq(bbox = 'Saint-Petersburg') %>%
add_osm_feature(key = "admin_level", value = '5') %>%
osmdata_sf()
# Но есть проблемы с кодировкой
iconv(admin_osm$nodes$tags, from="UTF-8", to="UTF-8")
iconv(admin_osm$nodes$tags, from="UTF-8", to="cp1251")admin_osm[["osm_multipolygons"]]
Объединяем слои и строим свою первую карту!
# Устанавливаем одинаковые проекции для всех объектов:
isr_border <- st_transform(isr_border, 4326)
isr_city = st_transform(isr_city, st_crs(isr_border))
isr_icecream = st_transform(isr_icecream, st_crs(isr_border))
# Делаем карту (st_geometry() - вывести чистую геометрию, без других колонок или значений в объекте)
{plot(isr_border %>% st_geometry())
plot(isr_city, col = 'red', add = T)
plot(isr_icecream, col = 'blue', add = T)}
Ура :)
Теперь сохраним данные (здесь есть свои подводные камни):
# write_sf(isr_border, "isr_border.shp", append = T, layer_options = "ENCODING=UTF-8")Загрузим в качестве проверки и убедимся, что всё работает
# double_isr_border <- st_read("isr_border.shp")
# plot(st_geometry(double_isr_border))
Бонус
Построим интерактивный график с помощью очень простой по своему синтаксису библиотеки mapview. Туториал по ней можно найти здесь. mapview предназначен для карт небольшого размера.
ЗАПУСКАТЬ ОСТОРОЖНО!!!
(на Windows возможны вылеты RStudio)
library(mapview)
mapview(isr_border) +
mapview(isr_city %>% st_geometry(), col.regions = 'red', legend = FALSE) +
mapview(isr_icecream %>% st_geometry(), col.regions = 'blue', legend = FALSE)
Open street map
naturalearthdata.com/downloads/
download.geofabrik.de/
Один из многочисленных обзоров челленджа
и много, много других ресурсов…
Pebesma E., Bivand R. Spatial Data Science. 2020.
Lovelace R., Nowosad J., Muenchow J. Geocomputation with R. 2021
Dorman M. Using R for Spatial Data Analysis. 2021.
Простое гугление даёт как минимум 10 книг о пространственном анализе в R.